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ABSTRACT 



Protostars grow in mass by accreting material through their discs, and this ac- 
cretion is initially their main source of luminosity. The resulting radiative feedback 
^ C| heats the environments of young protostars, and may thereby suppress further frag- 

[ mentation and star formation. There is growing evidence that the accretion of material 

JL \ onto protostars is episodic rather than continuous; most of it happens in short bursts 

that last up to a few hundred years, whereas the intervals between these outbursts 
of accretion could be thousands of years. We have developed a model to include the 
^ , effects of episodic accretion in simulations of star formation. Episodic accretion results 

in episodic radiative feedback, which heats and temporarily stabilises the disc, sup- 
pressing the growth of gravitational instabilities. However, once an outburst has been 
terminated, the luminosity of the protostar is low, and the disc cools rapidly. Provided 
that there is enough time between successive outbursts, the disc may become gravita- 
tionally unstable and fragment. The model suggests that episodic accretion may allow 
. disc fragmentation if (i) the time between successive outbursts is longer than the dy- 

■ namical timescale for the growth of gravitational instabilities (a few kyr) , and (ii) the 

quiescent accretion rate onto the protostar is sufficiently low (at most a few times 
10~ 7 Mq yr _1 ). We also find that after a few protostars form in the disc, their own 
episodic accretion events shorten the intervals between successive outbursts, and sup- 
press further fragmentation, thus limiting the number of objects forming in the disc. 
We conclude that episodic accretion moderates the effect of radiative feedback from 
young protostars on their environments, and, under certain conditions, allows the for- 
k> . mation of low-mass stars, brown dwarfs, and planetary-mass objects by fragmentation 

j_j ■ of protostellar discs. 



Key words: Stars: formation - Stars: low-mass, brown dwarfs - accretion, accretion 
disks - Methods: Numerical, Radiative transfer, Hydrodynamics 



1 INTRODUCTION 

Stars form in the dense cores of molecular clouds and grow in 
mass by accreting material from their surroundings. Due to 
the initial rotation and/or turbule nce of the parent clo uds 
discs form around protostars (e.g. iTerebev et al]|l984t ). A 
disc grows in mass by accreting material from the enve- 
lope. The material in the disc spirals inwards and onto 
the protostar, provided that there is an efficient mecha- 
nism to r edistribute angular momentum o utwards in the 



disc (e.g. Cartwright fc Stamatellos! l2010l : lAttwood et alj 



120091 : IWalch et al.l 120121 ). Otherwise the disc may become 
sufficiently massive and cool to fragment, forming mainly 
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low-mass stars and brown dwarfs (e. g. IWhitworth et all 
120071 : IStamatellos fc Whitworth! [2009a). Angular momen- 
tum can be redistributed by gravitational ins tabilities (GIs) 
or magneto-rotational instab i lities (MRI) (ILin fc Pringld 
1 19871 : iBalbus fc Hawlevl Il99ll . Il998l : lLodato fc Ricdl2004| ). 
During this process gravitational energy is transformed 
into radiation due to viscous dissipation in the disc and 
at the accretion shock around the protostar. This radia- 
tion heats the region around the protostar and may sup- 
press subsequent f r agmentation and further star formation 
l|Offner et al.ll2009l ; lBatdl2009bl ). 

Until recently computational studies of star formation 
in molecular clouds ignored t he effect of r adiative feedback 
from young protostars (e.g. Batd 12009a! ). This was done 
merely due to computational constraints but has no phys- 
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ical justification, since during the initial stages after the 
formation of a protostar, the accretion rate onto it can 
be rather high, and therefore its luminosity output can 
also be rather high (potentially up to 100s of Lq, even 
for a 1 Mq protostar). In the last few years, more de- 
tailed computational models including the effects of radia- 
tive feedbac k from protostars have been developed. It has 
been found (|Krumholzl|2006l; lBatell2009bl ; lOffner et alpOoj 
lUrban et alJl20ld : iKrumholz et alJl20ld : lOffner et alfaiof l 



that radiative feedback from newly formed protostars can 
suppress fragmentation and thus inhibit further star forma- 
tion. In particular, discs around young protostars are heated 
and stabilised and therefore discs may not contribute signifi- 
cantly to the production of low-mass stars and br own dwarfs 
1 120091 : lBat3l2009bh . [Bate! (|2009bh fc lKrumholzl 



Offncr et al.l 



201ll ) have concluded that radiative feedback plays a criti- 



cal role in regulating the stellar initial mass function. 

These studies have assumed that the accretion of mate- 
rial onto young protostars (and therefore their energy out- 
put) is continuous. However, there is growing observational 
evidence that accretio n is episodic, i.e. it happens in rela- 
tively short outbursts (lHerbidll977l; lDopitdll978l; iReipurthl 
ll989l ; lHartmann fe Kenvonll 19961 ; IGreene et al.ll200fih . 

FU Ori- type stars llHerbiel 1 19771 : IGreene et all l200Sl : 
IPenevaet al.ll2010l :lo" reen et alJl201ll ) provide the first piece 
of this evidence. These stars exhibit large increase in their 
brightness (~ 5 — 6 mag). The increase in brightness hap- 
pens within 1 — 10 yr and may last for tens to a few hun- 
dred years (lHartmann fc Kenvonll 19961 ; IGreene et al.1 1200II 
IPeneva et alJl2010h . The accretion rates during these events 
can be up to ~ 5 x 10 -4 Mq yr -1 ; this means that a large 
amount of material is delivered ont o the protostar d uring 
these outbursts (a few 10" 2 M , e.g. lBell fc Lm]ll994l ). 

Additional evidence for episodic accretion comes from 
the knots seen in many outflow s from young protostars 
(e.g. lDopitalll97i ; lReipurth|[l989l 'l. It is believed that out- 
flows parallel to the rotation axis of a protostar are driven 
by accretion of m aterial from its disc (|Pudritz et al.ll2007l ; 
IShang et al.l 120(371 ) . Therefore, periodically spaced knots in 
such jets are indicative of episodic accretion. 

The case for episodic accretion is further supported by 
the observed luminosities of young protostars. These lu- 
minosities are much lower than expected if protostars ac- 
crete steadily; this is referred to as the l uminosity problem 
ilKenvon et al]|l99Cl : lEvans et al.l 120091 ; lEnoch et all 120091 ; 



iDunham et al.l 



2010). For example, if one assumes that by 



the end of the Class phase (~ 10' yr), half the protostar's 
final mass has been accumulated, then for a final mass of 
1 Mq , the mean accretion rate onto the protostar must be 
~ 5 x 10~ 6 MQyr _1 , and the mean accretion luminosity 
must be ~ 25 Lq ; this is more than an order of magnitude 
large r than the observed average (~ 2 Lq; e.g. lEnoch et alj 
120091 ) . Part of the accreti on energy (i) may be u sed to power 
outflows and/or winds l|Pudritz et al.l 120071 ; IShang et al.l 
l2007i ;l Hansen et al. 1120121 ). fii) may be absorbed by the pro- 
tostar ( Hartmann et al. 19971 ; ICommercon et al.l 120111. or 



(iii) may be used to di ssociate and ionisize the accreting 
gas (|Tan fc McK cc 2004). However, these effects are not ex- 
pected to significantly reduce the protostellar luminosity. 
Recent observations have revealed that Class and Class I 
protostars have compara ble mean luminosities (I Evans et al.l 
120091 ; lEnoch et all [2009); this suggests that protostars ac- 



crete for longer periods than previously thought, thus 
the luminosity problem may be less pronounced. Never- 
theless, episodic accretion is still ne eded to explain the 
low luminosities of young protostars ( Dunham et al.l |2010| ; 
lOffner fc McKeell201ll ; [Dunham fc Vorobvovll2012i ). 

The physical cause of episodic accretion is unclear. It 
may be (i) due to gravitational interactions with companions 
or passing stars that perturb the disc and enhance accre- 
tion of material onto the protostar (|Bonnell fc Basticn| [l992l ; 
iForgan fc Ricell2010h. fii) due to thermal instabilities occur- 
ring i n discs (jrlartmann fc Kenvonlll985l ; iLin fc Papaloizo"ul 
1985; Bell & Lin 19 94]), or (ii i ) due to gravitationa l instabil- 
ities l|Vorobvov fc Basull2005l ; iMachida et al.ll201ll ). 

Another interesting possibility is that episodic accretion 
is caused by the complementary role of MRI an d GI in di scs 
l|Armitage et all l200ll ; IZhu et all 120071 . l2009al lbl l2010bl fah. 
The outer disc is sufficiently cool that it is susceptible to 
GIs, and these result in gravitational torques that transport 
angular momentum outwards, thereby allowing material to 
spiral inwards towards the central protostar. In contrast, the 
inner disc is too hot for GI, but, if it becomes hot enough, the 
gas is thermally ionised to a sufficient degree that it couples 
to the magnetic field and the MRI is activated (e.g. IClarkel 
2009; iRice fc Armitagell2009r ) ; the MRI then transports an- 
gular momentum outwards, allowing matter to spiral in fur- 
ther and onto the protostar. As a result matter accumulates 
in the inner disc region, and steadily heats up until it is 
sufficiently ionised to activate the MRI, whereupon the ac- 
cumulated matter spirals into the central protostar, giving 
rise to an outburst. Once most of the accumulated matter 
in the inner disc has been deposited onto the central pro- 
tostar, the temperature of the remaining matter falls, the 
MRI sh uts off, and the acc u mulat ion of matter resumes. 

In IStamatellos et all (1201 ll ) we incorporated the 
IZhu et all (1200981 ) model, in which MRI and GI work com- 
plementary, into radiative hydrodynamic simulations of star 
formation. We found that episodic accretion can limit the 
effect of radiative feedback from protostars, and therefore al- 
low disc fragmentation, and the formation of low-mass stars 
and brown dwarfs. Here, we expand this work and study the 
effect of different model parameters; we find that episodic 
accretion may allow fragmentation, but the process is self- 
limiting. 

In Section 2, we describe the computational method 
used to study star formation in collapsing molecular clouds, 
and in Section 3 we describe in detail the phenomenologi- 
cal model that we have developed to include the effects of 
episodic accretion in hydrodynamic simulations. In Section 
4, we present and compare the results of models (i) with 
no accretion luminosity, (ii) with continuous accretion lu- 
minosity, and (iii) with episodic accretion luminosity. We 
also examine the effect of the different parameters of the 
episodic accretion model. In Section 5 we discuss the con- 
ditions under which episodic accretion allows or suppresses 
fragmentation, and finally in Section 6 we summarize the 
results of this work. 



2 COMPUTATION METHOD 



We use the SPH code seren l|Hubber et all l2011al lbh to 
treat the gas hydro- and thermo-dynamics. The code in- 
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vokes an octal tree to compute gravity and find neigh- 
bours, multiple particle timesteps, and a Leapfrog integra- 
tion sc heme. The code uses time -dependent artificial vis- 
cosity (|Morris fc Monaghanlll997 | ) with pa rameters a = 1, 
/3 = 2a, and a Balsara switch ( Balsar all 199a ). so as to reduce 
artificial shear viscosity. 

The chemical and radiative processes that regulate 
the gas temperature ar e treate d with the meth od of 
IStamatellos et al.l (|2007bT ) (see also lForgan et al.ll2009l ). The 
net radiative heating rate for particle i is 



dui 



4^sb ft 



4 

BGR 



(1) 



The positive term on the right hand side represents heating 
by the background radiation field, and ensures that the gas 
and dust cannot cool radiatively below the background radi- 



ation temperature T B 



is the Stefan-Boltzmann con- 



stant, Ej is the mass-weighted mean column-density, and 
R R (pi,Ti) and K p (pi,Ti) ar e suitably adjusted Ross eland- 
and Planck-mean opacities (IStamatellos et al.ll2007bh . The 
method takes into account compressional heating, viscous 
heating, heating by the background radiation field, and ra- 
diative cooling/heating. The method has been extensively 
teste d l|Stamatellos et al.l l2007bl ; IStamatellos fc Whitworthl 
2008). Despite the fact that under certain circumstances the 
method gives a poor estim ate of the column densit y through 
which SPH particles cool dWilkins fc Clarkdl2012h. it works 
well within high density fragments ( Young et al. 20121 '). and 
it reprod uces the macroscopic results of ID and 3D sim- 
ulations jMasunaga fc Inutsu ka 2000; Boss & Bodcn heimerl 



19791; [Whitchousc & Bate 20oi lBolev et al.| [2006; Cai et all 



2008) 



The gas is assumed to be a mixture of hydro- 
gen (70%) and helium (30% ) , with an equatio n of state 
jBlack fc Bodenheimerl Il975l ; iBolev et al.l 120071 ) that ac- 
counts for the rotational and vibrational degrees of freedom 
of molecular hydrogen, and for the different chemical states 
of hydrogen and helium. The dust and gas opac ity is set to 
K,(p, T ) — ko p a T b , where kq, a,b are constants (|BeIl fc Linl 
1 1994 ) that depend on the physical processes contributing to 
the opacity, for example, ice mantle melting, the sublima- 
tion of dust, molecular a nd H~ contributions (for details see 
IStamatellos et al.ll2007bl ). 

Up to densities p ~ 10 -9 g cm -3 , the self-gravitating gas 
dynamics, energy equation, and associated radiation trans- 
port are treated explicitly. Wherever a gravitationally bound 
condensation with p > 10 -9 gem -3 is formed, it is presumed 
to be destined to form a protostar. Then, in ord er to avoid 
very small timesteps, it is replaced with a sink (Bate et al. 
Il995l ), i.e. a spherical region of radius 1 AU. The sink inter- 
acts with the rest of the computational domain only through 
its gravity and luminosity. Matter that subsequently flows 
into a sink, and is bound to it, is assimilated by the sink, 
and ultimately ends up in a protostar at the centre of the 
sink. Normally, a sink is considered to correspond to a star 
forming in the simulation. In this work, we assume that a 
sink corresponds to a star plus a part of the star's accretion 
disc, i.e. the disc inside the sink radius; this part of the disc 
will hereafter be called the "inner accretion disc" (IAD); it 
is the part of the disc that is not simulated with SPH, but 
is modelled phenomenologically (see Section [3] for details). 

Radiative feedback from protostars that form during 



the simulation is taken into account by invoking a pseudo- 
background radiation field with temperature T PGn (v) (see 
Eq. [TJ that is a function of positi on relative to all the pro- 
tostars present in the simulation (IStamatellos et all [2007b; 
IStamatellos fc Whitworthll2009bl ). 



w 4 +E{ii^ 



(2) 



where L n and r n , are the luminosity and position of the n* 
protostar. At all times, the luminosity of a protostar is given 
by 



Mn 

M 



L©+/ 



GM n M n 



(3) 



where M n is the mass of the protostar, R„ its radius, and 
M n is the accretion rate onto it. The first term on the right- 
hand side is the intrinsic luminosity of the protostar (due to 
contraction and nuclear reactions in its interior). The sec- 
ond term is the accretion luminosity. / is the fraction of the 
accretion energy that is radiated away at the photosphere 
of the protostar, r ather than being expended e.g. driving 
jets and/or winds (|Pudritz et all 120071 ; IShang et all 120071 ; 
lOffner et ai1l2009l ). In this work we set / = 0.75. We also 
assume R n = 3 R(n is the t ypical radius of a young pro- 
tostar (|Palla fc Stahlerlll993l ). During the initial stages of 
star formation the accretion luminosity dominates over the 
intrinsic luminosity of the protostar. 

The above equations result in a background temper- 
ature profile around each protostar that is variable with 
time and dependant upon the accretion rate onto the proto- 
star and the protostar's mass. The background temperature 
drops away from the protostar as r -2 , which is the expected 
radial temperature profile on an envelope around a young 
protostar. 



3 A TIME- DEPENDENT MODEL OF 
EPISODIC ACCRETION 

Simulations of self-gravitating hydrodynamics on the scale 
of molecular cloud cores (i.e. sizes from 10 4 to 10 5 AU and 
masses from 0.3 to 10 Mq) can achieve sufficient resolution 
to capture the formation of discs around young protostars, 
and the effects of GIs i n the outer re g ions o f such discs (e.g. 
Goodwin etHI l2004al ; iBatel l2009bl . |2012| ; lAttwood etahl 
2009). However, such simulations cannot capture what hap- 
pens in the inner disc region (within a few AU from the 
central star). In this region numerical (artificial) viscosity 
is larger tha n the viscosity du e to the effect of GIs (e.g. 
Clarke! |2009l; IZhu et all l2009al ; ICartwright fc StamatelQ 
20101 ; iForgan et al.lboill ). and hence matter s driven to- 
wards the central protostar at an artificially enhanced rate. 
Additionally, to avoid short timesteps and computationally 
prohibitive running times, sinks on the order of 1 AU are 
invoked. It is normally assumed that any matter that enters 
the sink flows instantaneously onto the central protostar. 
However, at this scale magnetic fields could be important 
in determining the rate at which matter is delivered to the 
protostar, and their effect needs to be included in the mod- 
els. 

To take into account the physics that is dominant 
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in the inner disc region, without increasing greatly the 
computational cost, we have incor porated the resu lts 
of est ablished semi-analytic models (|Zhu et alJ l2009al fbl 
l2010bllal) into hydrodyn amic simulations of star formation 
( Stamatellos et al 1 l201ll ). We assume that the mass accreted 
into a sink is deposited initially onto an inner accretion disc 
(IAD) inside the sink, where it piles up until it becomes hot 
enough that thermal ionisation couples the matter to the 
magnetic field. At this point the MRI is activated, trans- 
porting angular momentum outwards, and thereby allowing 
the matter accumulated in the IAD to spiral inwards and 
onto the central protostar. Hence the sink mass is divided 
between the central protostar (M*) and its IAD (M IAD ): 



AL 



(4) 



Matter is assumed to flow from the IAD onto the protostar 
at a rate 



(5) 



where M BRG is the quiescent accretion rate, i.e. when there 
is no outburst, and the second term is a much higher ac- 
cretion rate that obtains only when MRI acts to transport 
angular momentum in the IAD. The IAD is presumed to 
be a continuation of the much more extended accretion disc 

outside the sink, which is sim ulated explicitly. 

In detailed disc models (|Zhu et al.l l2009al lbl. l2010bl lah 
the matter in the IAD couples to the magnetic field, 
due to thermal ionisation, once the temperature reaches a 
threshold of T„ RI ~ 1400 K. Using an Q-parameterization 



jShakura fc Sunvaevlll973h fo r the e ffective viscosity deliv- 
ered by the MRI. IZhu et al.l (|2010bl ) estimate that the ac- 
cretion rate during an outburst is 

M MRI ~5xlO- 4 M yr- 1 (2!S£l) , (6) 

where a^-g, is the effecti ve Shakura-Sunyayev parameter 
llShakura fc Sunvaevlll973l ) for the viscosity provided by the 
MRI. IZhu et all (|2010bf ) also find that the duration of an 
outburst is 



At, 



0.25 kyr (^f) ' 



M, 



O.2M 

1/9 



2/3 



10" 5 MQyr- 1 



(7) 



and onto the IAD. We assume that the critical temperature 
for the MRI to be activated is reached when enough mass 
for an MRI-enabled outburst has been accumulated in the 
IAD, i.e. 



M 1AU > M UM 



^MRI At MRI 



Using Eqs. [6] and [7] we obtain 

\ 2 / 3 



M IAD >O.13M 







M, 



10" 5 M yr" 1 



1/9 



(8) 



(9) 



The detailed thermal and ionisation balance in the IAD 
is not modeled explicitly since this is not the main concern 
of this pap er. Based on observations and models of FU Ori- 
type stars (|Hartmann fc Kenvonll99d ; IZnu et al.ll2010bh . we 
assume that rapid accretion onto the protostar, at a rate 



given by 



Mmri — 1.58 



At 



is initiated as soon as M IAD exceeds M MR1 ; t Q is the time 
at which this occurs, and the outburst terminates at t Q + 
At um . The factor 1.58 



e/(e — 1) is to allow for all the 
mass of the IAD to be accreted onto the star within At,,,,,. 



The time taken to accumulate (or re-accumulate) M N 



is At . 



At, 



! 13kyr 



AL 
0.2M<7 



2/3 



Af, 



1O- 5 M yr- 



-8/9 



(11) 



This is effectively the interval between successive outbursts, 
and, as will be explained later in the paper, crucially it is 
much longer than the duration of the outburst (see Eq. [7]), 
and significantly longer than the growth time of GIs (a few 
kyr). 

In this formulation there are two free parameters. The 
first is a MRI , which controls the strength and duration of the 
outburst; increasing a MRI makes for a more intense, shorter 
outburst. The mass delivered in each outburst, and the du- 
ration of the interval between outbursts, are essentially inde- 
pendent of Q MRI , which is fortunate, in the sense that a MRI 
is rather uncertain. O bservations and simulations sugges t 
q mri = 0.01 to 0.4 (e.g. iKing et al.ll2007l ; llsella et all 120091 ). 
The second free parameter is the quiescent accretion rate 
onto the protostar, M BRG , i.e. the rate at which material is 
assumed to accrete onto the protostar when the MRI is not 
active. 



4 THE ROLE OF EPISODIC ACCRETION IN 
LOW-MASS STAR FORMATION 

To determine the role of episodic accretion in star 
formation, we simulate the collapse of a star form- 
ing cloud core and the subsequent genesis of proto- 
stars. The initial prop ert ies of the molecular cloud core 
i|Goodwin et al.l l2004al lbl l200rj ; lAttwood et~aH 120091 ') are 
chosen to match the observed properties of prestel- 
lar cores dWard- Thompson et all 1 1994 IXndre et al.l Il996l; 
Ward-Thompson et aljfl999l: lAndre et al.ll200ol; lAlves et all 



200 ll; iKirk et al 



Stamatellos et al 



.20051 ; IStamatellos fc Whitworth! 12003 
2004). The initial density profile is 



p(r) 



/^KERNEL 



(12) 



KERNEL 



= 3 x 10 18 g cm 3 is the central density, 



Here p 

and -R KERNEL = 5, 000 AU is the radius of the central region 
within which the density is approximately uniform; the free- 
fall time at the centre is ~ 40 kyr. The outer envelope of the 
core extends to R r , 



,,. E = 50,000 AU, so the total core mass 
5.4 Mq, and the mass initially inside R KERNBL 
, KJ ,,. <FL — 1.1 Mq. The gas is initially isothermal at 
T = 10 K, and hence the initial ratio of thermal to grav- 



is M, 
is M„ 



CORE 



itational energy is q t 



U, 



THERM 



/\U C 



= 0.3 . We 



impos e an initial r andom, divergence- free, turbulent velocity 
field (|Batdl2009bh . with power spectrum Pkdk oc k 4 dk, to 
match the scaling laws observed in m olecular clouds (|Larsonl 
Il98ll ; iBurkert fc Bodenheimerl |2000|), and amplitude such 



that 



J\U G 



0.3 jjiiina et al.ll 19991 ). 
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Table 1. Summary of seven simulations of collapsing turbulent cores; with no radiative feedback from protostars (run eaO), with 
continuous radiative feedback (run eal), and with episodic radiative feedback with various parameters (runs ea2-ea6). omri is the 
viscosity parameter due to MRI, M BRG is the quiescent accretion rate onto the protostar, Nf is the number of objects (stars and brown 
dwarfs) formed in the simulation, Mi their masses (at the end of the simulations, i.e. at 100 kyr), and ti their formation time. For the 
simulations with episodic accretion, and hence episodic radiative feedback, we also list the percentage of a star's life during which it is 
experiencing high MRI-driven accretion (tEA,i/t*,i)' Multiple entries of Mi, ti, and tEA,l/t*,i refer to different objects formed, in order 
of their formation time. 



run id 






N f 


Mi (M Q ) 


ti (kyr) 


*EA,i/**,i (%) 


eaO 






1 


0.49, 0.32, 0.005, 0.23 


78.7, 87.2, 88.8, 90.0 




eal 






1 


0.59 


78.7 




ea2 


0.1 


io- 7 


2 


0.63, 0.37 


78.7, 89.9 


4, 6 


ea3 


0.01 


io- 7 


1 


0.60 


78.7 


41 


ea4 


0.3 


io- 7 


3 


0.63, 0.26, 0.12 


78.7, 90.7, 91.0 


1, 2, 1 


ca.o 


0.1 


IO" 8 


3 


0.56, 0.35, 0.11 


78.7, 89.0, 90.9 


5, 2 


ea6 


0.1 


10~ 6 


1 


0.64 


78.7 


4 



The molecular cloud core is represented by 10 6 
SPH particles, so each SPH particle has mass m SPH ~ 
5 x 10" 6 M Q . The minimum resolvable mass is therefore 



3x10 M Q . This is much smaller 



than the theoretical minimum mass for star formation 
(~ 3 x 10~ 3 Mp), usually referred to as the opacity limit, 
e.g. IWhitworth fc Stamat cllos 2006|), and therefore all self- 
gravitating condensations are well resolved. 

We perform seven simulations, all with the same initial 
conditions, and differing only in their treatments of the lu- 
minosities of protostars (see Table [T] Figs. [TJHJ) . In the first 
simulation (eaO) there is no radiative feedback from young 
protostars, and in the second simulation (eal) the proto- 
stellar radiative feedback is continuous. In the remaining 
simulations (ea2-ea6) episodic accretion is invoked, which 
leads to episodic radiative feedback. Each run uses different 
values for the MRI viscosity parameter, a MRI , and/or the 
quiescent accretion rate onto the protostar, Af BRG . 

In all simulations the initial collapse leads to the for- 
mation of a protostar at t~ 79 kyr, and this acquires an 
accretion disc. The simulations only diverge after this junc- 
ture, as the disc growth and evolution depend on how the 
protostellar feedback is treated. The simulations are run un- 
til 100 kyr. 

It is not trivial to define a disc forming in a hy- 
drodynamical simulation and to calcu late its properties 
(|Machida et al.ll2010l ; IWalch et al.ll2012l) . Firstly, we calcu- 
late the angular momentum vector of the sink/protostar, 
and assume that the disc plane is perpendicular to it. Then 
we define concentric rings around the sink, we assume that 
the disc has thickness of z = 0.2r (note that since there 
is little mass outside O.lr from the disc midplane the esti- 
mated disc properties are not significantly affected by the 
choice of the disc thickness), and calculate the (azimuthal) 
average of the radial (v r ) and azimuthal (v^) velocity. We 
also estimate the Keplerian velocity at the radius of each 
ring, using Vk(ri) = (GM;/r;) 1,/2 , where r; is the radius of 
the ring, and Mi the mass inside this radius (including the 
mass of the sink). We calculate the disc radius (and after- 
wards the disc mass) in two ways: (i) we estimate the most 
distant ring for which v r (ri) < 2-KV^ri), i.e. matter can go 
around the sink at least once before it falls onto the sink 
(Fig. [3] red lines), and (ii) we estimate the most distant 



ring for which v$(ri) > 0.9 Vk(ri), i.e. the rotation velocity 
is nearly Keplerian (Fig. [3] blue lines). Both criteria give 
comparable results for the disc mass and radius, although 
the disc is generally somewhat more compact with the latter 
criterion (see Fig. [3J). 

The details of the simulations and their results in terms 
of star formation are summarised in Table [T] Henceforth, 
we will refer to the first object formed in a simulation as 
the primary protostar, whereas the rest, which are formed 
in the disc of the primary protostar, for simplicity will be 
referred to collectively as secondary protostars (even though 
some of them have masses below the hydrogen burning limit, 
~ 80 Mj, and therefore are brown dwarfs and/or planetary- 
mass objects). 



4.1 No radiative feedback 

In this run (eaO) mass accretes continuously onto the pro- 
tostar (see Figs. [1] [2] top panels), but this does not trans- 
late into radiative feedback, as the luminosity of the proto- 
star in Eq. [3] is artificially set to 0. The disc around the 
primary protostar grows in mass (up to ~ 0.2 Mq; see 
Fig. [5] top) and bec omes larger (radius out to ~ 100 AU) 
(cf. IZhu et all [20121 ) . It becomes gravitationally unstable 



enough (iMatzner & Levinl 20051; 


Whitworth & Stamatellosl 


20061: IWhitworth et all 120071; 


Stamatellos et al.l l2007al; 


Kratter et all 20081 IStamatellos & Whitworth 2009a! bl). 



The temperature at a distance of 100 AU from the primary 
protostar, a characteristic radius at which fragmentation is 
expected, remains very low around 10 K (Fig. |4ji, red line). 
The Toomre parameter Q at the same radius drops to near 
unity as the disc grows in mass, indicating that the disc 
becomes gravitationally unstable (Fig. |4]d, red line). Three 
secondary objects form in the disc (see Table[T|; two of them 
are low-mass stars and the third is a planetary-mass object. 
The planetary-mass object is qui ckly ejected from the disc, 
as see n in previous simulations (|Stamatellos fc Whitworth! 
2009a); this is critical for term inating its growth by ac- 
creting material from the disc l|Stamatellos fc Whitworth] 
2009b) . The results presented here are consistent with previ- 
ous studies that, in the absence of protostellar radiative feed- 
back, report disc fragmentation and low-mass star forma- 
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Figure 1. Evolution of the accretion disc around the primary protostar forming in a collapsing turbulent molecular cloud core. The 
colour encodes the logarithm of column density, in g cm -2 . The different rows present time sequences from 85 to 91kyr for different 
treatments of the protostellar radiative feedback. Top row (eaO), no radiative feedback: the disc around the primary protostar increases 
in mass, becomes gravitationally unstable, and fragments to form 2 low-mass stars and 1 planetary-mass object. Middle row (eal), 
continuous accretion and continuous radiative feedback: the disc grows in mass, but radiative feedback keeps it sufficiently hot so that 
it does not fragment. Bottom row (ea2), episodic accretion and episodic radiative feedback: the disc becomes gravitationally unstable 
(first column) but the growth of GIs is suppressed by heating due to an accretion burst (second column); however, after this burst the 
disc cools again, GIs quickly develop, and eventually the disc fragments; one low-mass star forms in the disc. 



tion (|Bate et al . 2002; Goodw irTet al.ll2004al ; lAttwood et all 
2009; Bate 2009a; Offncr ct al. 2009). 



4.2 Continuous radiative feedback 

In this run (eal; Fig. [T] middle panel) mass accretes con- 
tinuously onto the protostar, and gravitational energy is 
radiated away according to Eq. [3] Therefore the protostar 
continuously heats its environment. The luminosity of the 
protostar varies from 10 — 200 Lq (Fig. [2J middle panel). As 
in the previous case the disc around the primary protostar 
grows in mass (up to ~ 0.4 Mq; see Fig. middle) and ex- 
tends in radius (out to ~ 100 AU). The disc is continuously 
heated by the protostar and therefore it is stabilised against 
fragmentation, despite its large mass. The temperature at 
100 AU remains relatively high at around 50 K at all times 
(Fig. [4^,, blue line), and the corresponding Toomre param- 
eter is always larger than 2 (Fig. 2p, blue line), indicating 
a stable disc. There are occasions when the disc becomes 
marginally unstable and spiral arms start developing in the 



disc; this leads to more material driven inwards in the disc 
and eventually onto the protostar, and therefore more en- 
ergy is released that heats and stabilises the disc, dampening 
any spiral st ructure. The disc is therefor e in a self-regulating 

stable state (|Lodato fc Ricdl2004l2005h . 

These results agree with the studies of lBatd <|2009bD & 
lOffner et al. (2009), which show that when protostellar ra- 
diative feedback is included, disc fragmentation is generally 
suppressed and fewer low-mass stars and brown dwarfs form. 



4.3 Episodic radiative feedback 

In this run (ea2; a MM =0.1, M BRO = HT 7 M vr" 1 ; Pig-ffl 
bottom panel) mass accretes onto the protostar episodically 
using the model described in detail in Section The mass 
flows into the inner disc region and then it accretes onto the 
protostar only when the conditions are right for the MRI to 
operate. The accretion rate is low (10 -7 Mq yr -1 ) when the 
MRI is inactive, and high (~ 10~ 3 Mq yr" 1 ) when the MRI 
is active. Therefore, the release of energy from the protostar 
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Figure 2. The effect of radiative feedback: (a) the accretion 
rates, (b) the luminosities, and (c) the sink/star/IAD masses 
of the primary protostar formed in the simulation, (i) with no 
radiative feedback (run eaO; top panel), (ii) with continuous ra- 
diative feedback (run eal; middle panel), and (iii) with episodic 
radiative feedback (run ea2; bottom panel), with « MRI = 0.1, 
bug = 10~ 7 Mo yr -1 . In the first two simulations 'sink' and 
'star' are equivalent, whereas in the episodic accretion run a sink 
comprises the star and its inner accretion disc. In the run with 
continuous radiative feedback the accretion luminosity stays high 
at ~ 10 — 200 Lq, whereas in the run with episodic radiative 
feedback due to episodic accretion of material onto the protostar, 
the luminosity stays low (SI Lq), apart from during the episodic 
events in which the luminosity rises up to ~ 10 3 Lq. 




100 



100 



85 90 
time (kyr) 



100 



Figure 3. The estimated disc masses and radii calculated using 
(a) v r (ri) < 2itv < i > (ri) (red lines), and (b) Ju(rj) > 0.9 v k (ri) (blue 
lines), for the run (i) with no radiative feedback (eaO; top panel), 

(ii) with continuous radiative feedback (eal; middle panel), and 

(iii) with episodic radiative feedback (ea2; bottom panel). The 
black lines denote the total mass of the (primary) sink in each 
simulation (for runs eaO and eal this corresponds to the primary 
protostar; for run ea2 this corresponds to the protostar plus the 
IAD). Star signs in between graphs mark the formation time of 
primary and secondary protostars. In the run with episodic accre- 
tion (ea2) the bands in between graphs mark the intervals where 
episodic accretion occurs (black: episodic accretion onto the pri- 
mary star, red: episodic accretion onto a secondary protostar). 
The estimates of disc masses and radii are reliable only until sec- 
ondary stars form in the discs. Episodic accretion events disrupt 
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Figure 4. The azimuthally-averaged (a) temperature (top), and 
(b) Toomre parameter Q (bottom) at radius 100 AU from the 
primary protostar plotted against time, for the run with no ra- 
diative feedback (eaO; red), with continuous radiative feedback 
(eal; blue), and with episodic radiative feedback (ea2; black). 



is also episodic; for most of the time it is below 1 Lq, but 
during episodic accretion events it can get up to ~ 10 3 Lq 
(Fig. [2] bottom panel). However, the protostellar luminosity 
does not suppress fragmentation, since it it high only for 
relatively short periods of time. 

The disc grows in mass and radius as in the case with 
no radiative feedback (mass up to ~ 0.2 Mq and radius 
out to ~ 100 AU; see Fig. [3l bottom), and the disc be- 
comes marginally unstable. At this instant (t ~ 86 kyr) an 
episodic accretion event occurs (see Fig. [1] bottom panel, 
second snapshot) disrupts the disc structure, and seems to 
suppress the growth of gravitational instabilities. However, 
this is only temporary; after the episodic accretion event 
terminates, the disc quickly becomes gravitationally unsta- 
ble again. This time the disc remains cool for long enough to 
fragment and produce a low-mass H-burning star (see Fig. [T] 
bottom panel, last two snapshots). 

At the characteristic radius of 100 AU the azimuthally 
averaged temperature (Fig. black line) is very similar to 
the corresponding temperature of the run with no radiative 
feedback (~ 10 K), except during episodic accretion events, 
when the temperature may go up to more than 100 K. 
Once the outburst terminates the disc quickly cools down to 
~ 10 K. Assuming that the disc surface density at 100 AU 
is E ~ 20 g cm -2 , the temparature T ~ 100 K, and the cor- 
responding Rosseland-mean opacity k ~ 1 cm 2 g~ 1 , using 
Eq. [1] we find that the cooling rate is u c ~ 50 erg s _1 g _1 . 
Additionally, assuming that the disc mass is Mdi sc ~ 0.1 Mq 
we find that the total specific internal energy of the disc 
is u = (3/2)(fcT//xm H ) ~ 6 x 10 9 erg g" 1 . Therefore, the 
timescale of disc cooling is t c = u/u c ~ 4 yr, i.e. much 
shorter than the interval between outbursts. 

The Toomre parameter Q is also similar to the case with 
no radiative feedback except during the episodic accretion 
events (Fig.[4p, black line). The interval between successive 
episodic accretion events is on the order of a few thousand 
years (~ 6 kyr), which is adequate time for gravitational 
instabilities to grow and the disc to fragment. 



4.4 The effect of different parameters for the 
episodic accretion model 

The episodic accretion model (Section has 2 free param- 
eters; (i) the viscosity parameter a MRI , which controls the 
strength and duration of the outburst; (ii) the quiescent ac- 
cretion rate M BRG , i.e. the rate at which material is assumed 
to accrete onto the protostar when the MRI is not active. 
These are rather uncertain. Observation s and simulations 
of di s cs suggest a„„ = 0.01 to 0.4 (e.g . iBalbus fc Hawlevl 
ll99StlKing et alj|2007l : llsella et alj|200g ). whereas observa- 
tions of young stars suggest quiescent accretion rates around 
10" 11 - 10~ 6 M n yr" 1 (e.g. ICalvet et alj|2004 iNatta et all 



120041 : iMuzerolle et al.ll2005l ; iMohantv et al.ll2005l ). In the fol- 
lowing sections we examine the effect of these parameters on 
the growth of gravitational instabilities and disc fragmenta- 
tion. 



44.1 a h 



= 0.01,0.1,0.3 (M BRG = lO^Mgyr" 1 ) 



The viscosity parameter a MRI controls how fast matter ac- 
cretes onto the protostar when the MRI is active; the higher 
q mri is, the faster angular momentum can be transported 
outwards in the disc, and the faster matter accretes onto the 
protostar. Hence, as there is a specific amount of mass that 
is delivered within each episodic outburst (M Mm ; cf. Eq. [S]), 
the outburst duration, Ai MRI , is shorter for larger a MRI (cf. 
Eq. [7|). On the other hand, the outburst is stronger, i.e. the 
accretion rate onto the the primary protostar and the lumi- 
nosity during the outburst are larger (Eqs. trUTH)) . 

These are evident in Fig. [5] where the accretion rate 
and luminosity of the primary protostar are plotted for 3 
runs (ea3, ea2, ea4; top to bottom) with a Mm =0.01, 0.1, 0.3 
respectively. In the run with a MRI =0.01 (ea3) the durations 
of the episodic accretion events are longer, and disc fragmen- 
tation is suppressed. This is because the disc remains hot for 
longer periods of time, whereas the intervals in between out- 
bursts are rather short and gravitational instabilities cannot 
fully develop (see Fig. [5]). The growth of the disc around the 
primary protostar is very similar to the case with continuous 
radiative feedback. 

In the run with a MRI =0.3 (ea4) the outbursts are more 
prominent but much shorter in duration; there is enough 
time for gravitational instabilities to develop and the disc 
to fragment. In this case two secondary low-mass H-burning 
stars form in the disc (see Table [1]). 

We conclude that the strength of the outbursts is rather 
irrelevant regarding disc fragmentation. The crucial factor 
that determines whether the disc around the primary proto- 
star fragments is the interval between successive outbursts. 
Shorter outbursts and longer intervals between them favour 
disc fragmentation. 



44.2 



M BRG = 10- 



10- 



TO^Mgyr" 1 (a MRI =0.1) 



The quiescent accretion rate M BRG determines how much 
energy is released by the protostar when the MRI is inac- 
tive; the higher the M BRG , the higher the luminosity released 
by the primary protostar, the hotter the attendant disc be- 
comes. 

This is shown in the luminosity plots for the three runs 
of increasing M BRG (ea5, ea2, ea6) in Fig. [7] For M BRG = 
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Figure 5. The effect of increasing the viscosity parameter 
q mri ■ ( a ) * ne accretion rates, (b) the luminosities, and (c) the 
sink/star/IAD masses of the primary protostar formed in the 
simulation with episodic accretion and quiescent accretion of 
M BRG = lO- 7 M yr- 1 , (i) with a Mm = 0.01 (ea3; top), (ii) 
with « MRI = 0.1 (ea2; middle), and (hi) with ct MRI = 0.3 (ea3; 
bottom). The duration of the episodic accretion events is longer 
for smaller a MRI , and, as a result, the intervals between successive 
events are shorter, suppressing fragmentation of the disc around 
the primary protostar. 




Figure 6. The azimuthally-averaged (a) temperature (top), and 
(b) Toomre parameter Q (bottom) at radius 100 AU from the pri- 
mary protostar plotted against time, for runs with different ct MRI 
(red: a MRI = 0.01, run ea3; black: a MRI = 0.1, run ea2; blue: 
a MEI = 0.3, run ea4). In these runs M BRG = 10 — 7 Mq yr~ 1 . 



10 8 Mq yr 1 the luminosity of the protostar when the MRI 
is inactive is below 0.1 Lq, whereas in the run with Af BRG = 
10 _6 MQyr _1 the luminosity is 1 — 2 Lq. As a result the 
disc remains hotter in the second case. In Fig.|8]we plot the 
temperature at 100 AU from the primary protostar; it can 
be seen that, when the MRI is inactive, in the former run the 
temperature remains low at about 10 — 15 K, whereas in the 
latter run the temperature is relatively high at ~ 20 — 30 K. 
This temperature is high enough to stabilise the disc and 
suppress fragmentation. 

In the run with M BRG = 10" 8 M yr -1 (ea5) the disc 
fragmets and two low-mass H-burning stars form, in the run 
with M BRG — 10~ 7 MQyr -1 (ea2) the disc fragments but 
only one low-mass H-burning stars form, and finally in the 
run with M BRG = 10 -6 Mq yr -1 (ea6) the disc is stabilised 
by the released energy and it does not fragment. 

We conclude that high quiescent accretion rate M BRG 
makes disc fragmentation less favourable, with M BRG > 
10 _6 MQyr _1 being enough to totally suppress fragmenta- 
tion. 



5 THE IMPORTANCE OF EPISODIC 

ACCRETION FOR DISC FRAGMENTATION 

Radiative feedback from young protostars has the potential 
to heat and stabilise their disc s, suppressing fragmentation 
ijOffner et aLll2009l : lBate1l2009bl ). Episodic accretion provides 
a way to limit the effect of protostellar radiative feedback. 
In this model the accretion of matter onto the protostars 
occurs not continuously, as previously has been assumed, 
but in short episodic outbursts. Therefore energy is released 
in short batches rather than continuously. 

The crucial factors that determine whether fragmenta- 
tion is promoted by episodic accretion are (i) the interval 
between successive outbursts, and (ii) the energy released 
by the protostar when the outbursts are absent (i.e. when 
the MRI is inactive). Fragmentation will occur if the disc 
remains cool enough for long enough that GIs have time to 
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Figure 7. The effect of increasing M BRG : (a) the accretion 
rates, (b) the luminosities, and (c) the sink/star/IAD masses of 
the primary protostars formed in the simulations with episodic 
accretion and a MRI = 0.1, (i) with quiescent accretion rate 
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1 (ea5; top), (ii) with quiescent accretion 

©y r ~ 



1 (ea2; middle), and (i) with quies- 
1CT 6 Mm yr _1 (ea6; bottom). The 



1CT 8 M Q yr 

brg = iO -7 ^ 
cent accretion rate M BRG = iu — " Mq yr 
fragmentation of the disc around the primary protostar is not 
suppressed in the first two cases, with three and two secondary 
objects forming, respectively. However, a continuous, quiescent 
accretion of M BRG = 10~ e MQyr~ 1 is enough to raise the lu- 



minosity of the primary protostar to 
suppress fragmentation. 



1 Lq and completely 
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Figure 8. The azimuthally-averaged (a) temperature (top), and 
(b) Toomre parameter Q (bottom) at radius 100 AU from the 
primary protostar plotted against time, for runs with different 
quiescent accretion rates M BRG : (red: M„ — in _8 u„™-l 



run ea5; black: M„ 



BRG * 

10~ & M r . 



BRG — 

run ea2 



10 _s Mq yr" 
blue: M„„„ 



10 MQyr , run ea6). In these runs a N 



= 0.1. 



grow. Thus, provided that the intervals between successive 
outbursts are longer than the dynamical timescale (which 
at ~ 100 AU is ~ 1 — 2 kyr), and the quiescient accre- 
tion rate onto the protostar M BRG is at most a few times 
lO^Mgyr -1 , then fragmentation is possible, and indeed 
we find that it happens (see Table UJ. 

This is demonstrated in Fig. [9] where we plot the in- 
terval between successive outbursts (taking into account all 
young protostars in the simulation, i.e. the primary proto- 
star and the secondary protostars formed by disc fragmen- 
tation). Long intervals (~ 6 kyr) favour disc fragmentation 
provided there is still enough mass in the disc; when the in- 
tervals become shorter further fragmentation is suppressed 
(ea2, ea4, ea5). In the last run (ea6) fragmentation does not 
happen as M BRG = 10 _6 MQyr _1 , and the disc remains 
relatively hot (~ 20 — 30 K) and stable, despite the fact 
that the intervals between successive outbursts are quite 
long (^6 — 7 kyr). 

Episodic accretion may also act against fragmentation 
after a few secondary protostars have formed. These proto- 
stars also exhibit episodic outbursts of accretion, leading to 
outbursts of luminosity, so if they remain in the vicinity of 
the primary protostar they also heat and stabilise its disc. 
This is shown in Fig. 1101 in which we plot the luminosity of 
all the protostars formed in the simulation (for those sim- 
ulations where more than one protostar is formed, i.e. ea2, 
ea4, ea5). We see that after the second protostar has formed 
the intervals between outbursts become much shorter and 
fragmentation is therefore suppressed; this effect is further 
enhanced by the formation of a third protostar, which has 
its own outbursts, and therefore makes the intervals between 
outbursts (from all protostars) even shorter. 

We conclude that episodic accretion may initially pro- 
mote disc fragmentation and the formation of a couple 
of secondary low-mass stars (and/or brown dwarfs) in 
the disc, but suppresses subsequent fragmentation, limit- 
ing the number of objects that can form in the disc (cf. 
IStamatellos fc Whitworthll2009al ) 
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Figure 9. The intervals between successive episodic events (in- 
cluding events from the primary and from the secondary pro- 
tostars) for runs ea2 to ea6 (top to bottom) plotted against 
time measured from the formation of the primary protostar. The 
coloured stars correspond to the times of formation of each one 
of the secondary protostars, and the horizontal dotted line to an 
interval of 2 kyr, which is on the order of the dynamical timescale 
of fragmentation. 
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6 SUMMARY 

We have developed a phenomenological model to include the 
effects of episodic accretion (leading to episodic radiative 
feedback) in radiative hydrodynamic simulations of star for- 
m ation. The model i s based on the semi-analytical approach 
of lZhu et all (|2010bt) . and describes the basic aspects of ma- 
terial being episodically accreted onto a protostar and t he 
resulting episodic energy release (jstamatcll os et alll201ll ). 

The model has two free parameters, namely the vis- 
cosity parameter Q MRI which controls how fast matter flows 
onto the protostar when the MRI is active, and M gRG which 
represents the quiescent accretion rate onto the protostar, 
when the MRI is inactive. These parameters are poorly con- 
strained and they may vary from protostar to protostar in 
different environments. We have shown that if the accretion 
of material onto protostars is episodic rather than contin- 
uous, then radiative feedback from young protostars may 
not be as effectiv e in suppressing di s c fragmentat ion as pre- 
viously thought |Offner et alj|2009l ; lBatefl2009bl ). provided 
(i) that the interval between outbursts is longer than the 
dynamical timescale for fragmentation, i.e. a few thousand 
years, and (ii) the quiescent accretion rate onto the proto- 
star is less than a few lO _7 M0yr _1 . The time interval 
between successive outbursts depends on the details of the 
activation of the MRI, on how fast mass is delivered onto 
the inner accretion disc, and consequently on how fast mass 
infalls onto the disc, and on how fast it spirals into the in- 
ner disc region (see Eq. [TTJ. These physical processes are 
not well understood. 

Adopting typical values for the model parameters we 
find that disc fragmentation is initially promoted, and a few 
low-mass secondary objects form; however, further fragmen- 
tation is suppressed due to outbursts of energy from the sec- 
ondary objects formed, limiting the number of objects that 
form in the disc. 

Environmental factors may be important when consid- 
ering the effect of episodic accretion. For example, in a 
dense stellar environment perturbations from fly-bys may 
stren gthen GIs in the disc (|Thies et al boirjlForgan fc Rice! 
l2010l ) and deliver more mass to the inner regions triggering 
episodic accretion events. Therefore it may be expected that 
in dense stellar environments or in the centres of stellar clus- 
ters more outbursts would lead to stable protostellar discs, 
and consequently fewer low-mass stars and brown dwarfs 
being fo rmed. However, f or the moment there is no evidence 
for this (|Riaz et al.ll2012l ). 

We conclude that episodic accretion may provide an 
effective way to limit the effects of radiative feedback from 
newly formed protostars on their environments and facilitate 
disc fragmentation and the formation of low-mass hydrogen- 
burning stars, brown dwarfs, and planetary-mass objects. 



Figure 10. The luminosities of all protostars formed in the sim- 
ulations where multiple protostars form (top: ea2, middle: ea4, 
bottom: ea5). After secondary protostars form they also exhibit 
episodic accretion events and thus they increase the temperature 
in their immediate vicinity, i.e. the disc of the primary protostar. 
Therefore after a couple of secondary protostars have formed fur- 
ther fragmentation tends to be suppressed. This is because the 
intervals between successive luminosity bursts, i.e. when fragmen- 
tation is likclv. arc sienificantlv shortened (cf. Figure l9l. 
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